;
; Regrid CLM Nitrogen-Deposition file (ndep) to a new resolution. 
; Works for both climatology monthly ndep files as well as yearly dynamic files.
; Uses input environment variables to determine resolution and other
; parameters to describe set of files to use
; (input simulation year, output resolution, or file type [climatology of yearly dynamic]).
; Uses bld query to get input filenames. Also uses env variable CSMDATA for location
; of input files.
;
;  Erik Kluzek
;  Dec/14/2007
;  $Id: ndepregrid.ncl 28400 2011-05-16 05:46:46Z erik $
;  $HeadURL;
;
begin
  ; ===========================================================================================================
  ;
  ; IMPORTANT NOTE: EDIT THE FOLLOWING TO CUSTOMIZE
  ; Edit the following as needed to interpolate to a new resolution.
  ;
  ; Input and output resolution and sim_year for input data
  ;
  resin          = "1.9x2.5";              ; Input resolution:  1.9x2.5 (normally NOT changed)
  gridfile       = getenv("GRDFIL");       ; Get input grid file to use  for output resolution
  res            = getenv("RES");          ; Get output resolution from env variable
  rcp            = getenv("RCP");          ; Get input representative concen. pathway
  sim_year       = getenv("SIM_YR");       ; Get input simulation year from env variable
  sim_year_range = getenv("SIM_YR_RNG");   ; Get input simulation year range from env variable

  ; END OF SECTION TO EDIT
  ; ===========================================================================================================
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

  ; Defaults if env variable is NOT set

  if ( ismissing(gridfile) )then
     gridfile = "default";      ; default means use query and other env parameters to get filename
  end if
  if ( ismissing(res) )then
     if ( gridfile .ne. "default" )then
        print( "Must set output resolution description ('RES' env variable) when you supply a gridfile name (GRDFIL)" );
        exit
     end if
     res      = "48x96";      ; Output Resolution: 1.9x2.5, 48x96 etc.
  end if
  if ( ismissing(rcp) )then
     rcp      = "-999.9";     ; Input representative concentration pathway (2.6,4.5,6,8.5,-999.9)
  end if
  if ( ismissing(sim_year) )then
     sim_year = "2000";       ; Input/Output simulation year: 1850, 2000
  end if
  if ( ismissing(sim_year_range) )then
     sim_year_range = "constant";       ; Input/Output simulation year range: constant, 1850-2000, 1850-2100, or 2000-2100
  end if

  ; END OF SECTION TO EDIT
  ; ===========================================================================================================
  ;
  ; Use above to get filenames, and various other information needed such as dates
  ;
  if ( rcp .eq. "-999.9" ) then
     rcpmsg = "";
  else
     rcpmsg = "rcp"+rcp+"_";
  end if
  csmdata  = getenv("CSMDATA");
  clmroot  = getenv("CLM_ROOT");
  querynml = "bld/queryDefaultNamelist.pl -silent -justvalue -options sim_year="+sim_year;
  querynml = querynml+",sim_year_range="+sim_year_range+",bgc=cn,rcp="+rcp;
  if ( .not. ismissing(csmdata) )then
     querynml = querynml+" -csmdata "+csmdata;
  end if
  if ( ismissing(clmroot) )then
     querynml = "../../"+querynml;
  else
     querynml = clmroot+"/models/lnd/clm*/"+querynml;
  end if

  sdate     = systemfunc( "date +%y%m%d" );
  ldate     = systemfunc( "date" );
  print( "query string="+querynml )
  filename    = systemfunc( querynml+" -res "+resin+" -var stream_fldfilename_ndep -namelist ndepdyn_nml" );
  gridinfile  = systemfunc( querynml+" -res "+resin+" -var fatmgrid" );
  if ( gridfile .eq. "default" )then
     gridfile    = systemfunc( querynml+" -res "+res  +" -var fatmgrid" );
  end if
  if ( sim_year_range .eq. "constant" ) then
     filenameout = "fndep_clm_"+rcpmsg+"simyr"+sim_year+"_"+res+"_c"+sdate+".nc";
  else
     filenameout = "fndep_clm_"+rcpmsg+"simyr"+sim_year_range+"_"+res+"_c"+sdate+".nc";
  end if
  ;
  ; Open files
  ;
  print( "Interpolate from file: "+filename );
  if ( systemfunc("test -f "+filename+"; echo $?" ) .ne. 0 )then
     print( "Input ndep file does not exist or not found: "+filename );
     exit
  end if
  nc       = addfile( filename,    "r" );
  print( "Use gridin file:       "+gridinfile );
  if ( systemfunc("test -f "+gridinfile+"; echo $?" ) .ne. 0 )then
     print( "Input gridinfile does not exist or not found: "+gridinfile );
     exit
  end if
  ncgi     = addfile( gridinfile,  "r" );
  print( "Use grid file:         "+gridfile );
  if ( systemfunc("test -f "+gridfile+"; echo $?" ) .ne. 0 )then
     print( "Input gridfile does not exist or not found: "+gridfile );
     exit
  end if
  ncg      = addfile( gridfile,    "r" );
  print( "Output file:           "+filenameout );
  if ( systemfunc("test -f "+filenameout+"; echo $?" ) .eq. 0 )then
     system( "/bin/rm -f "+filenameout );
  end if
  nco      = addfile( filenameout, "c" );
  ;
  ; Define dimensions
  ;
  latgg = ncg->LATIXY(lsmlat|:,lsmlon|0);
  longg = ncg->LONGXY(lsmlat|0,lsmlon|:);
  nlat  = dimsizes( latgg );
  nlon  = dimsizes( longg );
  dimnames = (/ "time", "lat", "lon" /);
  if ( sim_year_range .eq. "constant" ) then
     ntime0 = 0;
     ntime  = 1;
  else
     ntime0 = 0;
     ntime  = dimsizes( nc->YEAR );
  end if
  dsizes   = (/ ntime-ntime0, nlat,  nlon /);
  is_unlim = (/ True, False, False /);
  filedimdef( nco, dimnames, dsizes, is_unlim );
  ;
  ; Define vars and add attributes from original file
  ;
  vars = getfilevarnames( nc );
  i = dimsizes(vars) - 1
  do while ( i .ge. 0 )
     dimlist = getfilevardims( nc, vars(i) )
     filevardef (    nco, vars(i), typeof(nc->$vars(i)$), dimlist );
     filevarattdef ( nco, vars(i), nc->$vars(i)$ );
     if ( isStrSubset( vars(i), "NDEP_" ) )then
        nco->$vars(i)$@_FillValue = -1.e36;
     end if
     delete( dimlist );
     i = i - 1
  end do
  fileattdef ( nco, nc );
  ;
  ; Add global attributes to output file
  ;
  nco@history            = ldate+": Regrid from "+resin+" resolution to "+res+" by ndepregrid.ncl";
  nco@source             = "Regrid from:"+filename;
  nco@gridfile           = gridfile;
  nco@ndepregridVersion  = "$HeadURL: https://svn-ccsm-models.cgd.ucar.edu/clm2/branch_tags/cesm1_0_rel_tags/cesm1_0_3_n04_clm4_0_32/models/lnd/clm/tools/ncl_scripts/ndepregrid.ncl $";
  nco@ndepregridRevision = "$Id: ndepregrid.ncl 28400 2011-05-16 05:46:46Z erik $";
  ;
  ; Copy coordinate variables over
  ;
  ; for yearly dynamic files...
  if ( sim_year_range .eq. "constant" ) then
     ntime0 = closest_val( toint(sim_year), (/nc->YEAR/) );
     if ( ntime0 .lt. 0 .or. ntime0 .gt. dimsizes(nc->YEAR) )then
        print( "ntime0 is out of range = "+ntime0+" sim_year is = "+toint(sim_year) );
        exit
     end if
     ntime  = ntime0 + 1;
  end if
  ntimef = ntime-1
  print( "Time range: "+ntime0+" to "+ntimef );
  nco->YEAR = (/nc->YEAR(ntime0:ntimef) /);
  nco->time = (/nc->time(ntime0:ntimef) /);
  print( "Years = " );
  print( nco->YEAR );
  nco->lon= doubletofloat( (/longg/) );
  nco->lat= doubletofloat( (/latgg/) );
  lat     = (/ nc->lat  /);
  lon     = (/ nc->lon  /);
  lono    = (/ nco->lon /);
  lato    = (/ nco->lat /);
  areai   = doubletofloat( (/ ncgi->AREA /) );
  sumarea = sum(areai);
  areai   = areai / sumarea;
  areao   = doubletofloat( (/ ncg->AREA /) );
  sumarea = sum(areao);
  areao   = areao / sumarea;
  print( "areai sum = "+sum(areai) );
  print( "areao sum = "+sum(areao) );
  if ( abs(sum(areai)-sum(areao)) .gt. 1.e-3 ) then
     print( "Input and output areas are significantly different" );
     print( "Area-in:  "+sum(areai) );
     print( "Area-out: "+sum(areao) );
     exit
  end if
  ;
  ; loop over variables
  ;
  if ( nlon .eq. 1 )then
     Cyclic = False
  else
     dx = lono(1) - lono(0); 
     ; The last longitude should equal 360-(dx-first_longitude)
     ; So for example if first_longitude=0,  last will be 360-dx
     ;                If first_longitude=dx, last will be 360
     expLast = 360.0 - (dx - lono(0));
     if ( abs(lono(nlon-1) - expLast) .lt. 1.e-4 )then
        Cyclic = True;
     else
        Cyclic = False
     end if
  end if

  print( "Cyclic grid: "+Cyclic );

  do i = 0, dimsizes( vars )-1

     vardimnames = getfilevardims( nc, vars(i) );
     ;
     ; If variable is not one of the dimensions -- regrid it and write to output file
     ;
     if ( (vars(i) .ne. "YEAR") .and. (vars(i) .ne. "month") .and. (vars(i) .ne. "lat") .and. (vars(i) .ne. "lon") .and. (vars(i) .ne. "LAT") .and. (vars(i) .ne. "LON") .and. (vars(i) .ne. "time") )then
       print( "Write variable: "+vars(i)+" to output file" );
       ;
       ; If time dimension
       ;
       if ( vardimnames(0) .eq. "time" )then
          nt = 0
          do t = ntime0, ntimef
             vart  = (/nc->$vars(i)$(t,:,:)/);
             if ( any(ismissing(vart)) )then
                print( "There are some missing values on input" );
                exit
             end if
             varto = linint2 ( lon,  lat,  vart, Cyclic, lono, lato, 0 );
             if ( all(ismissing(varto(nlat-1,:))) )then
                print( "Copy last latitude over" );
                varto(nlat-1,:) = varto(nlat-2,:); 
             end if
             if ( any(ismissing(varto)) )then
                print( "There are some missing values on output no="+num(ismissing(varto)) );
                exit
             end if
             nco->$vars(i)$(nt,:,:) = (/varto/);
             vart  = vart*areai;
             varto = varto*areao;
             print( "Year: "+nco->YEAR(nt)+" Input w-sum: "+sum(vart)+" Output: "+sum(varto) );
             nt = nt + 1
          end do
       ;
       ; without time dimension
       ;
       else
          var  = (/nc->$vars(i)$/);
          if ( any(ismissing(var)) )then
             print( "There are some missing values on input" );
             exit
          end if
          varo = linint2 ( lon, lat,  var, Cyclic, lono, lato, 0 );
          if ( all(ismissing(varo(nlat-1,:))) )then
             print( "Copy last latitude over" );
             varo(nlat-1,:) = varo(nlat-2,:); 
          end if
          if ( any(ismissing(varo)) )then
             print( "There are some missing values on output" );
             exit
          end if
          nco->$vars(i)$ = (/varo/);
          var  = var*areai;
          varo = varo*areao;
          print( "Input w-sum: "+sum(var)+" Output: "+sum(varo) );
       end if
     end if
     delete( vardimnames );

  end do
  if ( isvar("varto") )then
     delete(varto);
  end if
  if ( isvar("vart") )then
     delete(vart);
  end if
  delete(areao);
  delete(areai);
  if ( isvar("varo") )then
     delete(varo);
  end if
  if ( isvar("var") )then
     delete(var)
  end if

  print( "================================================================================================" );
  print( "Successfully created output ndep file: "+filenameout );

end
